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Abstract. The inner boundary of protoplanetary discs is structured by the dramatic opacity changes at the 
transition from the dust-containing to a dust-frcc zone. This paper explores the variety and limits of inner rim 
structures in passively heated dusty discs. For this study, we implemented detailed sublimation physics in a fast 
Monte Carlo radiative transfer code. We show that the inner rim in dusty discs is not an infinitely sharp wall 
but a diffuse region which may be narrow or wide. Furthermore, high surface densities and largo silicate grains as 
well as iron and corundum grains decrease the rim radius, from a 2.2AU radius for small silicates around a 47L0 
Herbig Ae star typically to 0.4AU and as close as 0.2AU. A passive disc with grain growth and a diverse dust 
composition must thus have a small inner rim radius. Finally, an analytical expression is presented for the rim 
location as a function of dust, disc and stellar properties. 
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1. Introduction 

The inner regions of protoplanetary discs are the birth- 
places of terrestrial planets. Understanding their nature 
and diversity is key to our picture of the formation of 
planetary systems. 

This paper uses Monte Carlo radiative transfer mod- 
elling and detailed sublimation physics to explore the in- 
ner rim structures of passive dust discs, and their merits 
in reconciling theory with observation. 

T Tauri and Herbig Ac/Be stars have prominent 
infrared excesses in their spectral energy distributions 
(SEDs), due to the presence of protoplanetary discs. The 
existence of such disc:s is known obscrvationally through 
statistics of photometric colours (Adams et al. 1987; 
Rydgren & Zak 1987; Calvet et al. 1992; Lada & Adams 
1992), spectral fitting (Kcnyon & Hartmann 1987), imag- 
ing (McCaughrean & O'dell 1996; Burrows et al. 1996), in- 
terferometry (Mannings fc Sargent 1997), and is inferred 
from studies of the Solar System and exoplanets. Based 
on broadband SEDs, discs are classified as either self- 
shadowed or flaring (Kenyon & Hartmann 1987; Meeus 
et al. 2001; DuUcmond & Dominik 2004) and their tem- 
perature structure can be approximated with two layers 
to explain solid-state silicate emission bands (Chiang & 
Goldreich 1997). In the inner regions of discs, dust sub- 
limates. The location where starlight reaches an optical 



depth of unity in the dust is referred to as the inner rim. 

The dust sublimation temperature is around 1500K, so 
the inner rim radiates in the near-infrared. 

In this paper, we use a new Monte Carlo radiative 
transfer code to study the effects of various physical pa- 
rameters on the inner rim structure of protoplanetary 
discs. Our motivation is outlined in Section 2. Dust subli- 
mation, backwarming and other theoretical aspects arc 
introduced in Section 3 and the results are presented 
in Section 4. A discussion and conclusions follow in 
Sections 5 and 6. 

2. Observational and physical motivation 

Herbig Ae/Be stars typically show excess near- infrared 
(NIR, A « 1.25 . . . 7/im) emission, corresponding to tem- 
peratures of ^ 1500K. Most fractional NIR excesses, 
fNiR = Lnir/L*, form a rising histogram from fNiR ~ 0.02 
to 0.22 and are well reproduced by models including dust 
sublimation in a hydrostatic disc (DuUemond et al. 2001; 
Isella & Natta 2005, hereafter DDNOl and IN05). Larger 
excesses, e.g. fNiR = 0.23 and 0.25 for AB Aur and WW 
Vul (Natta et al. 2001), and fwiR = 0.32 for HD 142527 
(Dominik et al. 2003), have proven impossible to repro- 
duce with existing models without introducing unrealistic 
assumptions about the rim scaleheight or dust sublimation 
temperature. 

Herbig Ae/Be and T Tauri stars have been shown to 

follow a Rrim L^^^ law for inner rim radii measured in 
the NIR regime (Monnier et al. 2005; Millan-Gabet et al. 
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Fig. 1. A schematic of the inner few AU of a disc. The 
colour map indicates the temperature, from < 300K (dark) 
to > 2000K {white). At the optically thin radius, Cbw = 1, 
dust can begin to condense. The classical rim, at Cbw — 4, is 
a maximally backwarmed dust wall. The actual rim, i.e. r = 1 
surface, is between these limits and preceded by an optically 
thin zone. 

2007), thought to coincide with the dust rim. Most Herbig 
Ae/Be rims in this sample have been shown to be scat- 
tered within expected hmits (Monnier et al. 2005) for a 
radius determined by dust subhmation, consistent with 
an optically thin inner hole. Some Herbig Be stars, with 
Ij-i, > 10'^Lq, conform to models where the inner hole is 
optically thick in the midplane and the rim radii are con- 
sequently smaller. Gas undergoing accretion onto the star 
has been assumed as the inner hole opacity source in these 
cases. 

The presence of a smooth, extended hot emission com- 
ponent interior to the dust sublimation radius has been 
deduced for several systems from interferometric obser- 
vations (Eisner et al. 2007; Tannirkulam et al. 2008b, a; 
Kraus et al. 2008; Isella et al. 2008). A hot emission com- 
ponent in the inner hole has recently gained further sup- 
port from spatially resolved spectroscopy (Najita et al. 
2009; Eisner et al. 2009). 

Meanwhile, the physico-chemical parameter space of 
dust rim models has been little explored in terms of vari- 
ations in global inner disc properties (surface density pro- 
file, composition) and dust type (astronomical silicate has 
been the standard opacity and homogeneous spheres the 
particle shape). The complex geometry and high optical 
depth of protoplanetary discs can now be studied with fast 
Monte Carlo radiative transfer codes (Min et al. 2009), 
providing further incentive to incorporate more detailed 
physics to clarify the limits of dusty rim models. 

3. Methods 

We consider models of passive discs, with temperature 
structures dominated by irradiation from the central star, 
and no accretional heating. (For a broader review of disc 
models, we refer the reader to DuUemond et al. 2007). 



The modern paradigm in passive disc models was devel- 
oped by Chiang & Goldreich (1997, hereafter CG97). In 
this approach, the disc consists of a cool midplane and a 
hot surface layer. Such models are highly successful in re- 
producing the general form and solid state bands of young 
star SEDs, but not the NIR excess, as emission from the 
dust sublimation region is not considered. 

Any realistic dust species will sublimate at or below a 
temperature of Tgubi ^ 1500K, corresponding to an emis- 
sion maximum around 2/im. In a disc, this implies a radius 
inside which no dust can exist. Terminating the CG97 den- 
sity distribution at such a radius results in a dust wall that 
is exposed to direct starlight, becomes hot and "puffs up" 
(Natta et al. 2001, and DDNOl). Backwarming, explained 
later, pushes the sublimation location to larger radii. The 
DDNOl models explain the NIR excess, but the emission 
is very inclination dependent, an effect contradicted by 
observations. 

A semi-analytical treatment including a gas density de- 
pendent sublimation temperature was applied successfully 
in explaining observed SEDs and deriving important prop- 
erties such as characteristic grain sizes (IN05 and Isella 
et al. 2006). The rim surface was shown to be curved rather 
than wall- like, due to the dependence of the sublimation 
temperature on gas density (IN05). Isella et al. (2006) suc- 
cessfully fitted the SEDs and interferometric rim radii of 
four out of five Herbig Ae stars with IN05 models. In most 
cases, the best fits required astronomical silicate grains of 
> 1.2/im size. However, the modelled rim radii were larger 
than those obtained by Eisner et al. (2004) by up to a fac- 
tor of three, and simultaneously fitting the SED and rim 
radius of AB Aur was unsuccessful. 

In a passive disc, the structure of the inner rim is deter- 
mined by the stellar parameters, the surface density, cool- 
ing efficiencies and sublimation temperatures of the dust 
species, and by backwarming. Because the dust tempera- 
ture and hydrostatic structure are also intimately coupled 
in a non-trivial geometry, finding static, stable solutions 
in terms of dust temperature structure is a difficult task, 
for which Monte Carlo radiative transfer is well suited. 

Monte Carlo radiative transfer has been used to extend 
the IN05 models to two grain sizes in thermal contact 
(Tannirkulam et al. 2007). It was shown that grain growth, 
specifically the presence of grains of different sizes and 
settling of large grains, strongly curves the inner rim. A 
comparison with these results is made in Section 5. 

3.1. The bulk gas and vapour densities 

The total amount of gas in a disc is referred to as the bulk 
gas in this paper. The bulk gas density determines the 
amount of dust present if a fixed gas to dust mass ratio, 
e.g. fgd — 100, is assumed. Dust grains always contribute 
some particles to the gas phase. These do not apprecia- 
bly change the gas mass but do maintain the gas-grain 
equihbrium, with gas phase particles of the dust species 
sticking onto a grain and balancing its evaporation. If all of 
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Fig. 2. Sublimation curves of the species used in this work. The 
vapour density of the species is used. At low vapour densities, 
corundum is the most refractory, whereas at higher values iron 
is relatively more refractory. At least one of the two species is 
always relatively more refractory than olivine. The expression 
for the curves is given in Eq. 1, and the fitting parameters 
in Table 3.2. The IN05 fitting formula has been abundance- 
corrected for this plot, otherwise it would lie ~ 50K below 
our curve for any olivine vapour density (see Appendix B). 
For astronomical silicate and forsterite, the olivine sublimation 
curve was used in modelling. 



a given dust species must be in the gas phase to maintain 
this equilibrium, we call its density the vapour density. In 
reality, many gas-phase species may combine to maintain 
the gas-grain equilibrium of any solid species. As at most a 
few species exist simultaneously in our models, we gener- 
ally neglect this fact, with some exceptions (see Sec. 5.2). 

3.2. Dust sublimation 

The temperature of a dust grain is determined by 
its wavelength-dependent opacity, and by the inten- 
sity and spectrum of the radiation field responsi- 
ble for the heating. The cooling efficiency, given by 
e = Kp(Tdust)/Kp(T*) = Cabs(Tdust)/Cabs(T*), is a 
ratio between the Planck mean opacity of the dust species 
at its own temperature to that at the stellar temperature, 
or for a single grain the ratio of absorption cross sections. 
In an identical stellar radiation field, a grain with a large c 
will achieve radiative equilibrium at a lower temperature 
than a grain with a small e, and generally epsilon increases 
with grain size. 

For a grain to be stable, the particle flux excaping 
from its surface must be balanced by the flux coming onto 
it from the gas phase. If the temperature gets too high, 
the grain will sublimate. The first rim models assumed 
that this would occur at Tgubi = const, e.g. DDNOl used 
1500K. Generally, the sublimation temperature of a dust 
species is a function of the partial pressure of its gas phase 
component, T^ubi = f(P)- The ideal gas equation of state 
can be used to replace the partial pressure P of the species 
with a function of its temperature and its vapour density. 



Pvapour- Loosely following Zemansky (1968), the vapour 
density at saturation pressure is given by the Clausius- 
Clapeyron equation: 



log (Pvapour) = B 



Lsubl 



(1) 



where A, B and C are thermodynamical quantities that 
can be fitted using laboratory measurements^. 

We take C = — 1 and obtain A and B from a pub- 
lished table of evaporation temperatures with correspond- 
ing bulk gas densities and abundances (Pollack et al. 1994, 
P94). The resulting parameters for Eq. 1, alongside those 
obtained from other sources, are summarized in Table 3.2. 
In fitting the P94 data, the partial pressure of each species 
was expressed through its abundance and molecular mass, 
and the bulk gas density. Parameters from other sources 
were converted (see Eq.2 of Lamy (1974) and Eq.lO of 
Cameron & Fegley (1982)). The table also gives Tsuh\ at 
Pvapour = 10~^^g/cm^ and references for each species, with 
boldface denoting the species used in this study. 

Sublimation curves of the main species used in this 
work are shown in Fig. 2. At vapour densities of 
10~^^ - ~^°g/cm^, olivine sublimates at around 1300K, 
iron at 1400K and corundum at 1500K. Corundum must 
have an abundance < 1% of the total dust mass and thus 
will necessarily have a much lower partial density than 
olivine or iron. 

3.3. Backwarming 

If a dust grain cannot cool into empty space in every di- 
rection, its temperature increases due to backwarming. 
Geometrically, the backwarming factor Cbw can be de- 
fined as the ratio of a full An solid angle to the solid angle 
subtended by the empty sky seen by the grain. For each 
photon emitted into the sky area obscured by nearby mat- 
ter, a similar photon returns, hence the term "backwarm- 
ing" . An estimate of the temperature of a dust grain on 
the rim surface is obtained from the equation 



47rR: 



2 

rim 



47r 

Cbw 



(2) 



where the energy absorbed by a grain at a distance 
im from a star of luminosity is equated with that 



^ Thcrmodynamically, A = AH/ In (10), whore AH is the 
enthalpy of vapourization, and C = (cq'/R — 1). The molar 
heat capacity has been broken into constant and temperature- 
dependent components, c'p = c[" + c"'. In addition, B = Bi -|- 
log (k//imp), where Bi is given by 



Bi 



In (10) {i 



/o crdT-J^c^dT 
T2R 



dT-hl 



and I is an integration constant. Sometimes, for example in 
Smithells (1967) as referenced by Lamy (1974), the full form 
of Eq. 1 must be invoked to provide a sufficiently accurate fit 

to experimental data, liowcver most often, as in this paper, the 
assumptions B ^ B(T) and C = — 1 are made. 
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Species 


Chemical 


Tsubl 


A 


B 


Reference for A and B 


Reference for opacity 




formula 


[K] 


(Eq. 1) 


(Eq. 1) 






Astron. silicate 


MgFeSi04 


1300 


28030. 


12.471 


Pollack et al. (1994, P94) 


Draino & Lee (1984) 


Olivine 


MgFeSi04 


1300 


28030. 


12.471 


P94 


Dorschner et al. (1995) 


Forsterite 


Mg2Si04 


1150 


26091. 


13.418 


Cameron & Fegley (1982, CF82) 


Sorvoin & Piriou (1973) 


Iron 


Fe 


1400 


21542. 


6.6715 


P94 


P94 


Iron* 


Fe 


1150 


20686. 


9.1134 


CF82 




Corundum 


AI2O3 


1500 


40720. 


18.479 


CF82 


Koike et al. (1995) 


Orthopyroxene* 


MgSiOs 




30478. 


14.898 


P94 




Quartz* 


SiOa 




26335. 


11.184 


Schick (1960); Lamy (1974) 




Water ice* 


H2O 


150 


2827.7 


7.7205 


P94 




Troilite* 


FeS 


680 


155.91 


-4.9516 


P94 (note: T^ubi = const) 





Table 1. A summary of the dust properties. (* - not used in the present work, given for completeness.) 
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Fig. 3. A schematic midplane dust temperature profile of an ef- 
ficiently backwarmed rim. In front of the rim, the optically thin 
equilibrium temperature of the dust grains is plotted. Once this 

drops below the sublimation temperature, dust can condense. 
Generally, the sublimation temperature decreases with radius 
because of the increasing scaleheight and decreasing surface 
density. Between the dust destruction and classical rim radii, 
there is a region in which dust may condense but need not be 
optically thick. Dust in this region is warm, but not close to the 
sublimation temperature. The dust temperature peaks locally 
at the optically thick rim location, and then falls rapidly to a 
low value inside the disc. See also Fig. 1. 

emitted by it. assuming a cooling efficiency e. Rigorously, 
Cbw is wavelength-dependent and its value can be ob- 
tained by solving a complex radiative transfer and dust 
sublimation and condensation problem, which is done in 
the Monte Carlo code we employ. However, as we show 
later, the description in Eq. 2, while simple, is useful in in- 
terpreting the results of our more rigorous models, as well 
as other work on the subject. Other things being constant, 
the rim location scales as Rrim cx; C^^. 

Cbw = 1 yields the temperature of a grain in an op- 
tically thin environment, and when TsuU is considered 
it gives the smallest possible distance from the star a 
passive grain may orbit for longer than its sublimation 
timescale. This limit has been used by e.g. Tuthill (2001) 



and Monnicr & Millan-Gabct (2002) to obtain lower limits 
on disc inner radii. It must be kept in mind that such lim- 
its are derived under assumptions regarding grain opac- 
ities, sublimation temperatures and geometry In assum- 
ing the rim was an instantaneously optically thick vertical 
wall, DDN2001 adopted Cbw = 4. This illustrates how far 
backwarming can move the rim. As Monnicr ct al. (2005) 
pointed out, the rim is actually between the two limits, 
depending on over what distance the dust becomes opti- 
cally thick. The concept of such an optically thin region 
in front of the rim is illustrated by Figures 1 and 3. 

On the rim, Cbw > 4 may also be true for special 
locations. One can imagine this for a dust grain sitting at 
the far end of a tunnel which extends into the rim (see 
also Appendix C). 

3.4. Monte Carlo radiative transfer with MCMax 

Radiative transfer through the complex dust geometries 
we wish to consider requires a flexible method in terms 

of density structures and dust composition. In addition, 
the coupled nature of the dust temperature and the hy- 
drostatic structure also calls for an iterative solution. The 
required flexibility and speed can be obtained by using 
Monte Carlo radiative transfer. 

We use here the Monte Carlo radiative transfer code 
MCMax (Min ct al. 2009). Because of an efficient imple- 
mentation of the diffusion approximation for regions with 
high optical depth, this code is especially suited for fast 
computations of geometries involving high optical depths. 

The radiative transfer is based on the Monte Carlo 
scheme of continuous reemission developed by Bjorkman 
& Wood (2001) and optically thin regions are treated with 
the method of Lucy (1999). 

To model the inner rim, dust sublimation was im- 
plemented into MCMax according to Eq. 1. The density 

structure of the dust is solved for by iterating the radia- 
tive transfer, dust sublimation and recondensation, and 
vertical hydrostatic equilibrium. After each iteration, the 
spatial grid is rebuilt in order to properly sample the new 
density structure. This ensures that in each iteration the 
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radial resolution around r « 1 is sufficient for At through 
any cell to be less than a fraction of unity. 

Condensation, if it can occur, must be treated care- 
fully. An upper limit is to transfer all the gas phase dust 
into grains, but if done too quickly this may overheat the 
cell, forcing sublimation and thus a loop. In the rim re- 
gion, between Cbw = 1 and 4, the change in solid dust 
mass in any cell is weighted so that the change is smaller 
the closer the cell is to the local Tgubi- The gas fraction 
of any dust species is constrained to be radially and verti- 
cally monotonous in this region. Furthermore, the optical 
depth added radially in front of the r = 1 surface cannot 
exceed 0.1. Some physical solutions might be suppressed 
by our constraints, but we expect a proper treatment of 
these to require time-dependent modelling. For a detailed 
discussion of the numerical aspects of implementing the 
sublimation physics we refer the reader to Appendix C. 

In reality, the sublimation of grains will cause the 
grains to lose mass, and thereby decrease in size. For prac- 
tical considerations, and because we want to systemati- 
cally study the effects of varying the grain size, we keep 
the size of the grains independent of the sublimation state, 
but increase or decrease their number density. 

3.5. Disc parameters 

3.5.1. General properties 

The disc model is parametrized as follows: a radial dust 

surface density distribution E(R) is described by two 
joined power laws, R~pi and R~p^, extending from an 
inner radius Rin out to Rout- The power laws are joined at 
the location Rp. The total diist mass, Mdust) gas to dust 
ratio, and a dust composition are specified. 

In general, we assume Rin = 0.03AU, Rout = 200 AU, 
Pi = 0.0; p2 = 1.5 and Rp = 4AU. A flat inner disc sur- 
face density makes the interpretation of our results more 
straightforward, especially for the optically thin zones. 

Most of our models assume a dust mass of Mdust ~ 
1O~''M0. For a parameter study of the inner disc surface 
density, this is varied from IO-'^Mq to W~'^Mq (lO-^M® 
to lO^Me)- 

The gas to dust ratio is 100 everywhere, and the gas 
and dust scaleheights are coupled. All our models use the 
star previously employed by IN05, with T* = lOOOOK, 

= 47Lq, and M^ = 2.5Mq. 

The IN05 comparison discs are parametrized as Rin = 

O.IAU, Rout = 200AU and pi = p2 = 1.5, with homoge- 
neous spheres (HS) of astronomical silicate as the opacity 
model. 

3.5.2. Dust types 

The opacities of the dust grains arc sensitive to their size 
and shape. For modelling the effects of the grain shape on 
the opacities, we use the statistical approach (Bohren & 
Huffman 1983; Min ct al. 2003). This allows to take the 
grain shape effects into account properly using computa- 



tionally favorable grain geometries. We apply the distri- 
bution of hollow spheres (DHS; see Min et al. 2005). 

For the composition of the grains previous studies 
have frequently used the so-called astronomical silicate 
(Draine &: Lee 1984). This species is most likely a compos- 
ite of different materials with different sublimation laws. 
Its composition is unknown and we are thus unable to 
derive a proper sublimation law for astronomical silicate. 
Therefore, in this paper we focus mainly on well character- 
ized materials which have been measured in the laboratory 
and are known to exist in meteorites. 

We begin by exploring discs with various grain sizes 
and surface densities of amorphous olivine. Then, models 
with two olivine grain sizes, and with dust types of dif- 
ferent transparency, refractivity and cooling effiency are 
presented. 

The presence of small, ^ O.l/xm grains in discs is re- 
vealed by fitting observed solid state emission bands. We 
use a model with such grains as a reference point. Large, 
> 10/im grains have a higher cooling efficiency than small 
grains, and as grain growth is known to occur in discs, they 
are also of great interest in modelling. To study the effects 
of large grains on the inner disc structure, we present mod- 
els with 10 or 100/im olivine grains added to 0.1/xm parti- 
cles. Discs of 1, 10, 100 and 1000/um olivine are presented 
in Sections 4 and 5. 

Surface density influences the inner rim in several 
ways: it sets the sublimation temperature, influences the 
backwarming efficiency and determines over what radial 
distance any given optical depth is reached. We vary the 
dust surface density in the inner disc from = 10"* to 
10^ as a parameter study. 

In modelling two grain types, we assume they are 
not in thermal contact, unless the contrary is explicitly 
stated. This is to allow more efficiently cooling or refrac- 
tory species to exist independently of less efficiently cool- 
ing or more volatile species. Thermal coupling will arise 
from the Monte Carlo radiative transfer in sufficiently 
dense regions, where all species will be in radiative thermal 
equilibrium. Unless stated otherwise, abundance ratios of 
9999/1, 90/10 and 10/90 are used in studying grain mix- 
tures. 

Above ~ lOOOK, amorphous silicates anneal to form 
crystalline silicates. The amorphous component is repre- 
sented by olivine in our models. It has been found that the 
annealed silicates are often iron-poor and magnesium-rich, 
e.g. forsterite and enstatite. Observations (van Boekel 
et al. 2004) and the short timescale of crystallization 
(Lenzuni et al. 1995) suggest that the crystalline silicate 
fraction on the inner rim surface approaches unity. 

We only present lO/im forsterite grains in this paper. 
Given the expected high abundance of crystalline particles 
in the inner disc, we also perform full scattering computa- 
tions for our olivine and forsterite mixtures. Using a 0.1%, 
10% or 90% forsterite abundance allows to follow the tran- 
sition from an amorphous- to crystalline-dominated rim, 
assuming the former is represented by olivine and the lat- 
ter by forsterite. 
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Iron grains, plausibly a sink of iron left over from 
forsterite formation, is another component we explore. 
Small iron grains cool more efficiently than small olivine 
grains, furthermore iron is significantly more refractory 
in the high vapour density regime, as seen in Fig. 2. 
Abundance ratios of 9999/1 and 90/10 are used in our 
olivine and iron disc models. 

Corundum has been suggested to exist in discs (e.g. 
Lenzuni et al. 1995). It is very refractory in the low 
vapour density regime, and forms a "highly refractory 
envelope" around olivine together with iron in Tgubi- 
Pvapour parameter space, as seen in Fig. 2. We use a trace 
olivine/corundum mass abundance of 9999/1 and com- 
pare it with a model where this ratio is 99/1. The latter 
is a convenient limit, approaching one obtained by taking 
a mixture with cosmic abundances and putting all alu- 
minium into corundum. 

Materials such as olivine do not exist in single molecule 
form and decompose during sublimation (Duschl et al. 
1996). Thus, particles transferred to the gas phase from 
one grain species may contribute to upholding the gas- 
grain equilibrium of another species (Dominik et al. 1993). 
We have attempted to take this into account by using the 
summed partial pressures of amorphous olivine and crys- 
talline forsterite to uphold the gas-grain equilibrium of 
both. Any pure iron released by the sublimation of olivine 
is not considered as a separate dust species in these mod- 
els. 

For comparison with previous studies we also present 
several models using astronomical silicate, in these cases 
the olivine sublimation law is applied. 

4. Results 

The results of our disc modelling are presented in this 
section. We first put forward single grain type models 
with realistic opacities. Then, dependencies on grain size 
and surface density are discussed. Finally, the effects of 
adding various amounts of different grain sizes or types are 
shown. A detailed comparison with the results of IN05 and 
Tannirkulam et al. (2007) may be found in Appendix B. 

The optical depths arc computed using the dust Planck 
mean opacity at the stellar temperature of lOOOOK. 

4.1. Olivine grains of various sizes 

For the main batch of models, we used opacities of iron- 
containing olivines, with sublimation properties from P94. 

4.1.1. The developing rim 

The first important result is that the inclusion of detailed 
condensation and sublimation physics leads to a rim loca- 
tion which is not entirely stable. This is demonstrated in 
Fig. 4, where we show the location and sharpness of the 
rim as a function of the computational iteration number. 

We begin the computation by placing the rim at the 
Cbw = 1 location (iteration 0), followed by two backwarm- 
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Fig. 4. Distances along the midplane, in AU, where r = 0.1, 1 
and 10 is reached, versus iteration number. The evolution of 
an Edust = 25g/cm^ inner disc with 10/im olivine over 10^ 
iterations reveals the oscillatory long-term behaviour of our 
models. Dust is built up over ~ 150 iterations, after which 
the sublimation temperature is exceeded and a retreat occurs, 
during which the optically thin zone is narrow. The fractional 
extent of the rim radius variation is significantly smaller for 
smaller grains. 



ing iterations where the location is dictated by an upper 
limit on Cbw The radiative transfer (temperature deter- 
mination across the disc), vertical structure of the disc, 
and dust sublimation and condensation are iterated on 50 
to 200 times. For each iteration, we plot a horizontal bar 
in Fig. 4, showing the radii where the optical depth ra- 
dially along the midplane ranges from 0.1 to 10, with a 
mark on the r = 1 location. Thus a short, point-like ex- 
tent of the bar indicates a sharp rim that jumps from low 
to high optical depth very fast, i.e. the classical structure 
as proposed by DDNOl. 

Over 1000 iterations, we see that both the rim location 
as well as its inward motion are unstable. Around iteration 
150, the sublimation temperature is exceeded in a small 
region in the rim, which becomes very sharp, with the 
T = 0.1 and 10 locations moving close together. The back- 
warming efficiency increases and the rim recedes. Around 
iteration 250, a new cycle begins with optically thin zone 
formation. 

Similar variability in the rim location is evident in 
models with multiple grain types. The variability in rim 

structure which we have demonstrated can, in some cases, 
considerably affect the fractional NIR excess, fNiR, making 
it difficult to associate illustrative values with a particular 

model. We are working on clarifying the issue. 

It must be emphasized that while we have just de- 
scribed the changes in rim structure as if time-dependent, 
our modelling has no physical time dependence. Rather, 
our results seem to indicate that there is no unique, static 
solution for the dust distribution in the rim, at least when 
a single grain type is used, and variations should be ex- 
pected in a proper time-dependent treatment. 
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4.1.2. Grain size and surface density 

As we saw in Section 3, both the grain size as well as the 
surface density influence the rim location. This is demon- 
strated in the flrst three panels of the summary diagram 
in Fig. 5, which show the rim location and sharpness for 
models with three different grain sizes and a range of sur- 
face densities. For a given total (solid and vapour) dust 
surface density of Ig/cm^ (models 0-2), the rim of 0.1/im 
grains is at 2.22AU, while 10/im grains put the rim at 
0.43AU and 100/xm grains at 0.29AU. A shift of the rim 
location towards the star is also evident for surface den- 
sity variations with a fixed grain size of 10/im (models 
3-8) and of lOO/im (models 9-13). We will return to this 
point in Section 5 and show how this can be understood 
analytically. 

Another feature of the single grain type models is the 

relative sharpness of the rims. It seems that if enough 
condensible material is available, the rim will take a sharp 
form, as assumed in the classical rim models. The low sur- 
face density computation in model 3 reveals the situation 
in discs with low amounts of condensible material avail- 
able. In model 3, the r = 0.1 to 10 region extends from 
0.67 to 1.88AU because the opacity is too small for the 
optical depth to rise quickly with radius even if all con- 
densible material is put in the solid phase. A similar effect 
is seen in models 9-13, where the same parameter study 
as with models 3-8 is carried out with 100/im grains. 

While in these models, the low surface density or very 
large grain size might seem artifical, effects similar to those 
seen here also emerge in multi-grain models with trace 
species that might be present in the very hot inner regions 
of discs. We will now turn our attention to these models. 

4.2. Two sizes of olivine grains 

Increasing the grain size generally allows cooling to gain 
in efficiency compared to heating. We will now present the 
effects of this in models where small and large grains are 
mixed. 

Models 14-17 of Fig. 5 contain two different grain sizes 

of the same material, in this case olivine. How the optically 
thin region is cooled by the large olivine grains, and then 
filled up as the abundance is increased, is demonstrated 
sequentially as temperature maps in Fig. 6. 

In the upper row of Fig. 6, the O.l/zm grains (left panel) 
have a rim at 2.22AU, the lO/um particles (right panel) at 
0.43AU. Differences in dust temperature throughout the 
inner disc are evident. 

Model 14, seen in the middle row, is a disc with a 
0.1/im to 10/im grain mass ratio of 9999/1. The trace 
abundance of efficiently cooling large grains is sufficient 
to give rise to an extended, cool region of low optical 
depth, which shifts the condensation location of 0.1/im 
grains from 2.22AU (model 0) to 1.68AU. The r = 10 lo- 
cation coincides closely with this point, indicating that a 
large amount of 0.1/im grains exists, pushing up the opac- 
ity. An optical depth of r = 0.1 is reached at 0.81AU, 



indicating the optically thin region covers 50% of the rim 
radius and around 75% of the surface inside that radius. 
The effect of pulling in the rim is more pronounced in 

model 15, seen in the bottom row of Fig. 6, where 10% 
of the dust mass is in 10/xm olivine grains and the rest 
in O.l/xm grains. The inner rim is now at 0.45AU. The 
large grains dominate the structure, which is close to that 
of the disc with only large grains (model 1). The smaller, 
0.1/um grains exist all the way to the new rim location, 
in a region where they could never condense without the 
presence of the larger particles. 

Model 16 shows that a fractional mass abundance of 
10^'' of 100/im grains has very little effect on the small 
grain rim, for their low abundance is insufficient to con- 
tribute to the optical depth. Model 17 demonstrates that 
a 10% abundance of 100/im grains is sufficient for them to 
dominate the rim location at a total dust surface density 
of Ig/cm^. 

4.3. Olivine and forsterite 

Forsterite is expected to be very abundant in the in- 
ner rim region. We now present modelled disc structures 
obtained with no scattering (our usual assumption) and 
with isotropic scattering for mixtures of lO/xm olivine and 
forsterite. 

Since relative to olivine, forsterite grains are highly 
transparent to stellar radiation, one might imagine them 
to exist very close to the star. In addition to having a low 
heating efficiency in the stellar radiation field, however, 
forsterite grains cool very inefficiently at temperatures of 
1300. . .1500K. This causes them to sublimate at approx- 
imately the same radius as the iron-rich silicates, leading 
to no difference between the rim structure of model 1 and 
models 18 to 20 in the case with no scattering, as seen 
in Fig. 5. Including isotropic scattering (lower, dark gray 
lines of models 18 to 20) slightly increases the rim radius 
because of the increased backwarming efficiency. Note that 
at temperatures of a few hundred kelvin, the forsterite 
grains will cool much more efficiently and be much cooler 
than the iron-containing silicate grains. 

4.4. Small amounts of effienctly cooling or 
ultra-refractory material 

Small, 0.1/im grains of iron and corundum are more ef- 
ficiently cooling than olivine of the same size, and may 
exist significantly closer to the star. The high refractivity 
of corundum and iron provides another mechanism that 
allows grains to survive near the star. 

A mass abundance of 0.0001 in 0.1/im iron grains pulls 
the rim of similarly sized olivine from 2.22AU to 1.63AU 
and gives rise to an optically thin zone covering w 10% 
of the rim radius (model 21). While 0.1/im iron is cooler 
than olivine, it does not cool the inner disc as efficiently as 
corundum does (see below and Fig. 7). A 10% abundance 
of small iron pulls the rim to 0.71AU, but now the rim is 
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Fig. 5. r = 0.1, 1 and 10 radii of selected models. Total dust surface densities are given, and normalized abundances are 

tabulated. Cases with no scattering arc given {black, upper lines), with isotropic scattering indicated for forstcrito (18-20, dark 
gray, lower lines). The horizontal bars indicate the t = 0.1 and 10 radii. Models ... 2 provide a comparison to surface density 
(3 . . . 13), size (14 . . . 17) and composition (18 . . . 26) variations. 



very sharp (model 22), and the structure resembles that 
obtained with a 90/10 mass ratio of small and large olivine 
grains (model 15, see Fig. 6) 

Adding an abundance as large as 10% of either 0.1 
or 10/um iron particles to lO/xm olivine (models 23 and 
24) has little effect on the rim location (model 1), as the 
cooling efficiency of larger olivine grains is similar to that 
of iron. 

Model 25 shows the effects of a mass abundance of 

0. 0001 of 0.1/xm corundum grains in a disc of olivine of 
the same size. The corundum is significantly cooler than 
the olivine, as seen in Fig. 7. The rim has moved from 
2.22AU (the 0.1/um ohvine rim in model 0) to 1.57AU, 
with an optically thin zone beginning already at 0.84AU, 

1. e. covering « 45% of the radial distance and 70% of the 
surface area inside the r = 1 radius. 

Allocating corundum a mass fraction of 0.01 (model 
26), a rough upper limit based on the assumption that 
the cosmic abundance of aluminium is all in corundum, 
demonstrates that this species can draw the rim in con- 
siderably. Small, O.ljLtm corundum grains draw the rim 



in more efficiently than iron of the same size, both be- 
cause they are relatively more refractory at the low densi- 
ties they have in these models (as low as 10~^^g/cm^, see 
Fig. 2), and because they cool more efficiently at temper- 
atures typical of dust sublimation. 



A c;oniparison of how small abundances of O.l/ini 
corundum and iron in a disc of identically sized olivine 
cool the inner disc, and how their refr activity allows them 
to exist closer to the star, is presented as temperature 
maps of models 25, 21 and 22 in Fig. 7. Corundum (top 
row) cools more efficiently than iron (e.g. middle row), 
and with a further contribution from its high refractivity, 
it moves the rim to 1.57AU, closer to the star than an 
identical mass abundance of iron, which is at 1.63AU. A 
10% mass abundance of small iron grains is sufficient for 
them to dominate the rim structure and decrease the ra- 
dius considerably, to 0.71AU. Due to the similar cooling 
properties of 0.1/im iron and lO/xm olivine grains, models 
15 and 22 are similar. 



Kama, Min, Dominik: The inner rim structures of protoplanetary discs 



9 



< 



N 



< 



N 



< 



N 








2 3 

R [AU] 



4 







2 3 

R [AU] 



4 



5 



Fig. 6. Temperature (main) and density (insets) maps of disc inner regions for various olivine mixtures. The axes show radial 
distance from the star along the midplane and vertical height from it. The colour levels represent the dust temperature in steps 
of lOOK, for reference we give solid lines at 500K and 1300K. Also shown are the radial (dashed) and vertical (dotted) t — 1 
surfaces for the stellar radiation. The inset maps have total solid dust density contours at factors 2.7, 10, 10^, 10^ and 10^" below 
the maximum. The dust surface density in these models is Eg « Ig/cm^. A typical midplane dust density is 3 ■ 10~^'^g/cm^. In 
models with several dust species, the species are thermally decoupled, but radiative thermal equilibrium arises naturally deep 
inside the disc. Non-uniformities in the temperature maps reflect Monte Carlo noise and the resolution of the angular grid. Top 
row: Models of discs with one grain size. Shown is the dust temperature for discs with 0.1 (left panel) and lO^m (right panel) 
olivine grains (models and 1 in Fig. 5). Middle row: A disc with a small to large grain mass ratio of 9999/1 (model 14). 
Given are the temperatures of the small, 0.1/im grains (left) and the larger, 10/^m grains (right). Bottom row: A disc with a 
mass abundance ratio of 90/10 in small/large grains (model 15). We again give the temperatures of the O.l/im grains (left) and 
10/im grains (right) that co-exist in this disc. 



5. Discussion 

5.1. Modelling summary 

We have successfully employed the Monte Carlo radiative 
transfer code MCMax (Min et al. 2009) in modelling the 
inner rim structures of dusty passive protoplanetary discs. 

Results in Section 4 show that, especially for larger 
grains, a rigorous treatment of dust sublimation and con- 
densation, together with other rim physics in MCMax, 
may not always lead to a stable solution for a solid dust 
density and temperature distribution in a passive, static 
disc framework. 



5.2. Where is the rim? 

The rims of MCMax disc models are generally between the 
optically thin destruction and the backwarmed wall radii 
(the Cbw = 1 and 4 radii, respectively). Grains can begin 
to attenuate the stellar radiation field from the Cbw = 1 
location outward, allowing the rim to exist closer than a 
fully backwarmed estimate. That the rim location could 
vary like this has been pointed out earlier (e.g. Monnier 
et al. 2005), but we have explored the phenomenon in 
detail for the first time. 

For a given dust species, larger grains will generally 
cool more efficiently and thus will set the rim location, but 
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Fig. 7. Temperature (main) and density (insets) maps of disc inner regions for various grain mixtures. See Fig. 6 for a detailed 
legend. Top row: A disc with an olivine to corundum mass abundance ratio of 9999/1. Shown are the temperatures of the 
0.1/im olivine (left) and the 0.1/im corundum (right) grains. (Model 25 of Fig. 5, compare with model 0.) Middle panel: A 
disc with an olivine/iron mass abundance ratio of 9999/1 (model 21). Shown are the temperatures of the O.l/xm olivine (left) 
and iron (right) grains. Lower panel: A disc with an olivine/iron mass abundance ratio of 90/10 (model 22). Given are the 
temperatures of the 0.1pm olivine (left) and iron (right) grains. 



only if they are present in sufficient abundance, as shown 
in Figures 5 and 6. Large grains of oHvine can exist much 
closer to the star than small grains. A small, 0.0001 abun- 
dance of 10/im grains is enough to create an extended, 
cool optically thin region in front of the rim, which is still 
determined by 0.1/im grains. Allocating a fraction 0.1 of 
the mass to large grains is sufficient for them to determine 
the rim location. 

It was noted by IN05 that the most refractory species 
would determine the rim location, if it was also able to 
build up enough optical depth. They also pointed out that 
for one species, the largest grains would position the rim 
under the same assumption. We add the generalization 
that refractivity and cooling efficiency combine to deter- 
mine the species which sets the rim radius. A very high 
cooling efficiency may allow a species which does not have 
the highest sublimation temperature to determine the rim 
location. 



For a given cooling efficiency, the species with the high- 
est sublimation temperature will determine the rim loca- 
tion, again assuming a sufficient abundance of it is present. 
Corundum has been considered an ultra-refractory con- 
densate in discs, which it is if one assumes that all species 
have equal partial pressures. However, if the dust vapour 
densities are taken equal, as in vertical slices in Fig. 2, 
corundum is more refractory at vapour densities Pvapour < 
10~^°, but at Pvapour > 10"^*^, the sublimation tempera- 
ture of iron is higher and increases rapidly. Hence, iron 
is the most likely dust species to be responsible for rim 
temperatures above 1500K ^. 

^ This is subject to the validity of the thermodynamic data 
for iron from P94. Sublimation temperatures more than 200K 
lower than this are found with data from Cameron & Fegley 
(1982). However, P94 considered a protoplanetary nebula envi- 
ronment, where iron-rich silicates and other species contribute 
to the vapour pressure of iron. They assumed a silicate to iron 
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Fig. 8. Midplanc rim locations and dust surface densities for 
various sizes of olivine. Boxes represent the final computed 
MCMax iteration, usually between 50 and 200. The left solid 
line connects models with 100/im grains, and the right 0.1/um. 
Dash- dotted lines represent analytical computations for the 
same grain sizes with Eq. 3 for Cbw = 1 and a dust temper- 
ature of 1500K. The dotted line connects models with 10/im 
grains, and the dashed one 1000/im. Grain size variations: 
The rim radius Rrim changes by almost an order of magnitude 
for olivine grain size variations. Surface density variations: 
For a given grain size, increasing the surface density decreases 
Rrim, in one limit as a power law (shown by the long-dashed 
line). Models for a given grain size follow irregular lines due to 
variations between iterations (see Fig. 4). 



The above points to the fact that if grain growth pro- 
ceeds steadily and the dust has a complex composition fea- 
turing efficiently cooling and refractory metal oxides such 
as iron or corundum, the rim will exist close to the star, at 
0.4AU or closer for a 47Lq Hcrbig star if > 10"g/cm^ 
(see Figures 5 and 8). Evidence for grain growth is pre- 
sented in e.g. van Bockcl ct al. (2004) and Hcrbst et al. 
(2008). Relatively large rim radii, for a star of the type 
used here Rrim > 2AU, are likely to arise as the inner 
disc is cleared of gas by the magneto-rotational instability 
(Chiang & Murray-Clay 2007) or a planet, lowering the 
sublimation temperature. Dust may also be dynamically 
cleared. 

Surface density controls the rim location for a given 



grain size. This is described by a Rrim cx: S' 



-1/24.4 



dust power 
law in the sublimation-controlled regime and is illustrated 
in Fig. 8. The low density ends of lines in Fig. 8 branch off 
from the power law as there is insufficient mass to reach 
T = 1 at the dust evaporation location. An analytical de- 
scription of the rim location as a function of surface den- 
sity for both the sublimation-controlled and optical depth 
controlled regimes is given in Eq. 3. It is a curve which 
branches from the power law and asymptotically becomes 
parallel to the x-axis. 



mass ratio similar to ours, therefore we use their results with 
confidence but state that a full treatment of the multi-species 
gas-grain equilibria is desired. 



c. 



r^rim — 
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e 



cm 



(3) 

This equation, derived in Appendix A, yields the 

curves seen in Fig. 8. k! is the dust opacity at 0.55/im, 
in cm^/g. Note that the surface density profile is an- 
chored to Esubi, its value at the sublimation radius, with 

a radial dependence of E cx R~p. The properties of the 
central star are hidden in Si = M" "^" ^-0.252 .^-^^ 

_ l + 0.078p l+8.080p 

S2 = ^ * . The constants are given by 

Ci = 4.1982 • 10-6 and C2 = 4.84 • 10^. 

The good correspondence of the analytical curves de- 
scribed by the above equation with our numerical results 
shows that Eq. 3 captures the processes that are thought 
to dominate the rim location in dusty discs, and we pro- 
pose it as an aid in interpreting more sophisticated nu- 
merical models, but also observations. 

The sublimation properties of olivine were used in the 
above derivation. To use the presented formalism for any 
dust species, a power law fit to Pollack et al. (1994) data 
should be obtained (Tgubi = Gk (Pgas / [Ig/cm^])''', see 
also IN05 and compare with Eq. 1). 

Omittances in the above analytical approach include 
the temperature dependence of the cooling efficiency, 
e(Tdust), and the fact that partial condensation of dust 
occurs in front of the r = 1 location, i.e. the rim is always 
somewhat diffuse. These two effects are the main culprits 
in the discrepancies between the numerical and analytical 
curves seen in Fig. 8. 

We use a flat dust surface density profile E oc R-p, 
where p = 0, in the inner disc, however Eq. 3 is generic 
and allows to use the canonical p = 1 . . . 1.5. Furthermore, 
the sublimation location for a given surface density and 
grain type is independent of the surface density power law. 

Used with observed rim radii, Eq. 3 may prove useful 
in estimating the surface density in the inner disc. For 
this, the equation should be solved numerically for Sgubh 
inserting the stellar properties and adopting constraints 
on the grain opacity and the surface density power law. 

5.3. Is the rim diffuse or sliarp? 

The optically thin inner disc region is an important con- 
cept. A r < 1 region between Cbw = 1 and 4 is supported 
by any of three conditions: 1) a low surface density, 2) 
a high surface density but very low opacity, and 3) early 
stages of dust condensation in a cool region. 

Condition 1 can be met by a low-abundance of a 
species which is more refractory than the bulk of the dust, 
such as corundum (see model 25 in Fig. 5), or which can 
cool relatively efficiently, such as large grains or iron (mod- 
els 14 and 21 in Fig. 5). 
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Condition 2 could be met by a dust type which is very 
transparent in the optical, but cools efficiently in the NIR. 
If the scattering phase function of such a species is not 
strongly forward-peaked, a population of grains closer to 
the star than the rim (determined in this case by the sub- 
limation of another, more opaque species) could further 
give rise to an extended hot radiation zone through scat- 
tering. 

Condition 3 is met as a transient phenomenon in our 

models, which allow dust to slowly condense into an ini- 
tially empty inner disc, creating and then filling in an 
optically thin region, as illustrated by Fig. 4. Transient 
heating events such as a powerful flares may destroy dust 
in an inner disc region, which will subsequently go through 
a similar cycle of a transient extended optically thin zone 
during dust re-condensation. It would be surprising to ob- 
serve long-lived sharp rims composed of only large grains. 
However, our modelling also indicates that adding several 
types of dust to a model makes the rim structure much 
more stable. 

5.4. Dust differentiation 

Differentiation of dust types can occur in the optically thin 
region, where highly refractory or very efficiently cooling 
species can exist independently of more volatile or less 
efficiently cooling species. Maintaining such an extensive 
optically thin region in a disc with a diverse and broad dis- 
tribution of grain sizes and compositions requires consid- 
erable fine-tuning, thus sublimation-based differentiation 
in a static disc is unlikely. 

5.5. Observational implications 

The smallest observed characteristic rim radii for L^, k, 
50Lq stars are around 0.2AU (Millan-Gabet et al. 2007). 
These come from ring model fits to visibility curves. In the 
framework of our models, such small radii require grains 
of around lOOjivsi to be present at dust surface densities of 
Ed « 10°- -^g/cm^ (Fig. 5, in particular models 11 to 13). 

It may be feasible to put a lower limit on the inner 
rim surface density Sj-im of a system with Eq. 3, under the 
assumptions that a measured Rrim corresponds to Cbw ~ 
1, and that Riim is determined by the sublimation of grains 
with e « 1. The former is suggested by comparison of 
the Monte Carlo and analytical results in Fig. 8, and the 
latter reflects the assumption that large, efficiently cooling 
grains will dominate the inner disc. One could then use a 
measured Rrim with the stellar luminosity and mass to 
obtain a simple lower limit on the surface density. 

Furthermore, it is interesting to assume that the 
largest grains needed to satisfactorily fit the mid- to far- 
infrared SED of a system provide a lower limit on the size 
of the grains setting the location of the inner rim. Under 
this assumption, an upper limit for Erim could be obtained 
by using the cooling efficiency e of these grains in Eq. 3. 
If we further assume that these grains dominate the rim 



location instead of e « 1 grains, we are again left with a 
lower limit on Erim- 

Comparing these limits with observed disc masses and 

canonical E oc power laws will constrain the 

difference between the inner and outer disc surface density 
profiles. 

Our results hold if the optical depth of gas in the inner 
hole to stellar radiation is Tg,* <C 1 and, relatedly, if accre- 
tion can be neglected. If the gas continuum opacity in the 
inner hole is Tgas,* ~ 1, shielding will decrease the dust 
destruction radius by [exp(— 1)]^/^ Ki 20%, if the gas it- 
self does not approach stellar temperatures. A hot gaseous 
component may give rise to NIR emission observed inside 
the dust destruction radius in some systems (Eisner et al. 
2007; Tannirkulam et al. 2008b,a; Kraus et al. 2008; Isella 
et al. 2008). 

As Najita et al. (2009) have pointed out, modelling in- 
dicates that gas in the inner hole of a Herbig Ae/Be star 
should produce spectral features of CO and II2O, but in- 
stead a continuum process was found necessary to explain 
the compact NIR excess of MWC480. Potentially relevant 
opacity sources inside the classical dust sublimation ra- 
dius include H2O, H~, free-free emission and highly refrac- 
tory or transparent dust species (e.g. Najita et al. 2009). 
However, even the most refractory species hypothesized to 
exist in discs, such as corundum, are not expected to have 
a significant solid fraction inside the classical silicate dust 
sublimation radius. 

If a dust species is responsible for a extended hot emis- 
sion component reaching from a few stellar radii to the 
optically thick dust rim, it is likely to have low absorp- 
tion in the visible range and a high cooling efficiency in 
the NIR. Moreover, species with small absorption cross- 
sections and relatively isotropic scattering phase functions 
in the visible or NIR regimes may play a role in creating 
an extended radiating zone in the inner disc by scattering 
stellar or rim emission. 

Forsterite meets some of the above requirements, how- 
ever it is not able to cool efficiently in the NIR and thus 
differs little from olivine in terms of rim structure, as seen 
in Fig. 5. 

5.6. Looking ahead 

To make full use of our work, observed SEDs and interfer- 
ometric visibilities (as opposed to rim radii derived from 
them under various assumptions) need to be fitted simul- 
taneously with MCMax or a similarly capable code. A 
paper focussing on the observational applications of our 
models is in preparation. We are performing a study of 
the extreme values of observables that can be obtained 
with static dust models, as well as simultaneously fitting 
the SEDs and visibilities of specific objects. 

We are actively exploring whether the fractional NIR 
excesses of our models provide a source of additional infor- 
mation about the inner disc surface density and/or grain 
type. As fNiR varies with iteration as well as grain, disc 
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and stellar properties, a better handle on the variability 
of the rim is needed before this fraction is useful. 

Rim temperatures obtained from NIR SED-fitting and 

limits on inner disc gas densities would increase the use- 
fulness of Eq. 3 in constraining the amount and type of 
dust in the rim region. 

Wc support the conclusion, overlapping in part with 
Monnier et al. (2005), that detailed modelling of multi- 
species condensation and sublimation processes, as well as 
simultaneous treatment of radiative transfer in the gas and 
dust components, is needed to make further progress on 
the numerical side, and encourage gas and dust modellers 
to join efforts. 



6. Conclusions 

We have demonstrated a wide range of possible inner rim 
structures in a parameter study of grain size, composition 
and dust surface density, as well as multiple grain type 
models. To do this, we implemented dust sublimation and 
condensation physics into a fast Monte Carlo radiative 
transfer code. Here, we outline our main conclusions: 

1. The inner rim in dusty discs is not an infinitely sharp 
wall. Backwarming effects combine with sublimation 
and condensation in leading to a diffuse region that can 
extend over a significant fraction of the rim radius. In 
our treatment, dust generally begins to condense near 
the Cbw = 1 location and the r = 1 surface is closer to 
the star than in previous models. 

2. The inner disc surface density of a given species is an 
important parameter, because it determines the high- 
est possible temperature where the species can be sta- 
ble. High surface densities move the rim closer to the 
star. 

3. The dust component (species or size) which determines 
the rim location by building up an optical depth of 
unity closest to the star need not have the highest 
sublimation temperature, it may instead cool very ef- 
ficiently. 

4. If large particles are abundant enough to produce opti- 
cal depth, then the rim location will be dominated by 
silicates. If large partic;les are not abundant enough, 
corundum and iron grains will likely set the rim loca- 
tion. In either case, for our standard star, the rim lo- 
cation is typically around 0.4 AU, much smaller than 
models for 0.1/xm silicate grains would predict. 

5. For any given grain material and size, the rim location 
can be analytically expressed function of stellar 
properties and of surface density. We give this expres- 
sion in full form for a silicate-dominated rim. 

6. The optically thin region can cover 70% of the surface 
inside the rim radius in cases where a very efficiently 
cooling or highly refractory species is present in low 
abundance compared to the less efficiently cooling or 
more volatile bulk of the dust. 



7. In the case of very low surface densities, the build-up 
of optical depth will be slow, and the rim region will 
be very diffuse. 

Appendix A: Analytical relation between rim 
radius and surface density 

The rim location for a given grain type is a function of 
surface density at the rim location, as seen in Fig. 8. At 
high surface densities, this dependence is a power law. 
This power law relation is a limiting case of a more gen- 
eral relation describing the r = 1 location as a function 
of surface density. At low surface densities, the rim loca- 
tion moves away from the dust sublimation radius as not 
enough dust is available there to reach an optical depth of 
unity. 

The calculations that follow use the cgs system of 
units. This is important to keep in mind when using the 
final expressions, where everything is implicit except the 
relevant characteristics of the dust, the disc and the star. 

A.l. Sublimation-controlled regime 

We assume initially that the rim is at the dust destruction 
location, Rrim = R-subii E^nd wish to describe this radius as 
a function of dust surface density, Ed — Sgas/fgd, where 
fgd = 100 is the gas to dust ratio. Assuming the minimal 
rim radius is at the dust sublimation location, with Cbw = 
1, we first express the temperature of a grain at the rim 
radius from Eq. 2: 

where Cbw is the backwarming factor, e the c;ooling 
efficiency of the grain and the luminosity of the star. 
This expression assumes the disc is optically thin until the 
rim radius. 

Further assuming that the gas and dust are in thermal 
equilibrium, their scaleheights will be the same and we can 
relate the midplane dust density pd to the dust surface 
density via 

Sd(Rrim)=y^ Pd(Rrim)exp (^-^^dz, (A.2) 

where z is the vertical distance from the midplane and 
h is the scaleheight, h oc T^^R^j^^^. Integrating Eq. A.2 
and replacing the temperature with Eq. A.l, we obtain a 
relation connecting Sd(Rrim)5 Pd(R'rim) and Rrim- 

Using the IN05 power law fit to the sublimation law of 
olivine, we can write 

/ { \ 0-0195 

Td = T,„bi = 2000 3 - (A.3) 

Vad • lg/cm3 J 

Here, the bulk gas density, used by IN05 in fitting, 
has been broken into the dust density p^, abundance of 
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the species aa and the gas to dust ratio fgd- Taking the 
appropriate abundance from P94, we obtain the relation: 



R-subl — Ci 



J]0.04^/f0.02 



Cbw 

e 



0.53 



cm 



(A.4) 



This power law, where Ci = 4.198-10"^, describes the 
minimum radius where dust can condense as a function 
of the surface density at that location. Note that if one 
assumes full condensation of the dust, backwarming will 
move the rim further, an effect which can be explored 
by varying Cbw A more general expression for the rim 
location is formulated below. 

A.2. Optical depth controlled regime 

The inner rim is defined as the radial r = 1 location for 

a stellar photon at A = 0.55/im. For high surface densi- 
ties and under the assumptions of full condensation and 
Cbw = 1) this corresponds closely to the dust destruc- 
tion radius. However, the lower the surface density, the 
longer the distance along the midplane that photons have 
to travel to reach an optical depth of unity. Thus, at rela- 
tively low surface densities, one predicts (and our numeri- 
cal results show) a turn-off from the power law of Eq. A.4. 

We begin by integrating the optical depth through dust 
radially along the midplane: 



/3d(R)dR. 



(A.5) 



R.,ub 



Here, Rgubi is the dust destruction radius of Eq. A.4, 
Rrim is the rim radius, for which one assumes r = 1, and 
k' is the mass absorption coefficient of the dust at 0.55/im. 
Pd(R.)) the radial run of dust density, is obtained by mul- 
tiplying Eq. A. 2 with a factor Ssubi(Rsubi/R)~''- This is 
a surface density power law fixed at the sublimation lo- 
cation, i.e. Sgubi and Rgubi obey Eq. A.4. Assuming some 
power p relates the rim radius to the surface density at 
the sublimation location. 

Employing Eq. A.4 for Rsubi and inserting physical 
constants, we are left with the following description of 
the rim location: 



c. 



l + 4p 



C 



bw 



-0.252 



5.-l-0.078pg 



Cbw 

e 



l+S.OSOp s -1 



cm 



(A.6) 

This equation yields the curves seen in Fig. 8. P re- 
lates to the surface density power law as P = (1 -I- 4p)/4 
and the properties of the central star are hidden in 

l + 0.078p l+8-080p 

Si = M2-"i" L-0-252 ^nd S2 = M;^^ LT^. 
The constants are given by Ci = 4.198 • 10~^ and 
C2 = 4.84 • 10^. It is important to keep in mind that the 
sublimation properties of olivine were used in the above 
derivation. 
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Fig. B.l. Radial t — 1 results for the IN05 disc models {dashed 
lines) and the same discs modelled with MCMa^x {solid lines). 
Grain size goes from left to right as 1.3 and 0.1/im. A system- 
atic difference in radius is evident. 



Appendix B: Differences from IN05 and 
Tannirkulam et al. (2007) 

In modelling discs parametrized identically with those of 
Isella & Natta (2005, IN05), we find systematically smaller 
rim radii. The main reasons arc that the IN05 formalism 
yields a higher optically thin dust temperature and a lower 
sublimation temperature than the treatment used in this 
work, as discussed below. A similar study by Tannirkulam 
et al. (2007, TOT), which we also comment on, found re- 
sults different from both this work and IN05. 

We find rim radii smaller than those of both IN05 and 
TOT and for all modelled grains, i.e. 0.1, 0.5 and 1.3/Ltm 
astronomical silicate. As seen in Fig. B.l, the relative dif- 
ference increases with grain size, from AR^im ~ 10%Rino5 
for 0.1/im grains to ARrim ~ 35%Rino5 for 1.3/im. 

The IN05 disc model uses the radiative transfer equa- 
tion given by Calvet et al. (1991) for a semi-infinite slab to 
get a radial temperature structure inside the rim, and then 
calculates the vertical structure as in CG9T and DDNOl. 
A gas to dust ratio of 100 is used, and the total disc mass 
is distributed according to a E oc R~^-^ surface density 
power law extending from Rjn = O.IAU to Rout = 200AU. 
The dust opacity determines the temperature and verti- 
cal structures, and the gas density is used to compute the 
dust sublimation temperature and thus governs dust sub- 
limation. 

To obtain a law for the dependence of the dust sub- 
limation temperature on gas density, IN05 fitted Eq. A. 3 
to data from P94. The fitted points were Tgubi of olivine 
(olivine is used interchangeably with astronomical silicate 
in this appendix, as the two are assumed to differ in opti- 
cal, and not in sublimation properties) and Pgas, the bulk 
gas density. Because the grain-gas equilibrium is main- 
tained by the partial pressure of olivine, this approach 
yields a Tgubi = T(pgas) law which is applicable only to 
gas of the same composition as that of P94. By applying 
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this law to a disc with a purely olivine dust composition, 
IN05 obtained lower sublimation temperatures than us, as 
discussed next. 

P94 supplied the fractional mass abundances of the 
species in their gas. For olivine, it is 0.00264, i.e. 26.4% 
of the dust mass, assuming a gas to dust ratio of 100. 
Thus, in the IN05 approach, olivine contributes 100% of 
the dust mass, but the sublimation temperature is com- 
puted from the bulk gas density under the assumption 
that there is 100/26.4 times less olivine present, leading 
to a correspondingly lower Tg^bi. This underlies our choice 
to reduce the sublimation law to a dependence of T^ubi on 
the vapour density of the species, Pvapour, and not the 
bulk gas - by doing so, we may use the law for all nebular 
compositions^. 

Fitting the P94 data for the olivine gas density re- 
moved its abundance from our T^uU ^s^w, allowing us to 
compute the olivine sublimation temperature for any mass 
abundance. In an identical disc with pure olivine dust and 
fgd = 100, we obtain higher sublimation temperatures 
than IN05, because for the same bulk gas density, pure 
olivine can maintain equilibrium at a higher temperature 
than that given by the law IN05 used. 

The above is illustrated by Fig. B.2, which gives vari- 
ous temperatures in the rim region of an IN05 comparison 
disc with 0.1/Ltm grains, modelled with MCMax. The dash- 
dotted line is the sublimation temperature computed and 
used in MCMax, according to our fit to P94 data. The 
dotted line, seen to lie ~ 50K below the previous, shows 
Tsubi calculated from the same density structure using the 
IN05 sublimation law. 

To use the IN05 law generally, a "dummy" bulk gas 

density can be calculated from an olivine density using 
both the P94 olivine abundance of 0.00264 and fgd = 100. 
This density can be used in Eq. A. 3 to obtain Tgubi. We do 
this when using Eq. A. 3 in Appendix A, obtaining results 
which match our numerical models very well. 

For radiative transfer, IN05 adopted the approach of 
Calvet et al. (1991). The disc is assumed to be an irradi- 
ated semi-infinite slab. A temperature structure derived 
from this overestimates the disc temperature, because a 
semi-infinite slab c;an only cool in one direction, but a 
disc has a finite height. Two dust temperature profiles 
are given in Fig. B.2, one from MCMax (solid line) and 
the other computed for the same density structure us- 
ing Eq.l of IN05 (dashed line). The IN05 approach over- 
estimates the temperature inside the rim compared to a 
Monte Carlo computation. The backwarming factor as de- 
fined in Sec. 3.3 is useful in understanding this, as we will 
now discuss. 



^ This is vaUd under the assumption that olivine is the only 
species contributing to its own gas pressure. While this is not 
the case as there is no gas phase olivine molecule, it is the 

best we can do at present if we wish to model discs with a 
composition different from the P94 nebula. 
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Fig. B.2. Midplane dust temperature structures for an 
MCMax model of the IN05 disc with O.l/im astronomical sil- 
icate grains. The MCMax dust temperature structure (solid 
line) is seen to drop more rapidly, and is significantly cooler, in 
the rim region than the estimate obtained for the same density 
structure with the radiative transfer equation used by IN05 
(dEished line). The ~ 50K difference between our estimate of 
the dust sublimation temperature (dotted line) and that of 
IN05 (dash-dotted line) is also clear. In the MCMax model, the 
rim radius is at 1.05AU. The intersection of the extrapolated 
optically thin temperature curves and the T^ubi computed by 
IN05 is at ~ 1.16AU. Inside the rim, the MCMax dust tem- 
perature and that computed using the IN05 radiative transfer 
equation coincide better in later iterations, when more dust 
has condensed there and made the rim structure more similar 
to a wall. 

The expression derived by IN05 (their Eq.6) from C91 
for the optically thin dust temperature is 

T\r, =0)^Tt= (2f, + n, (B.l) 

where Tq is the dust and Tgtar the stellar temperature, 
r the distance from the star, R^, the stellar radius, // is the 
sine of the incidence angle of starlight on the rim, and e 
the cooling factor as defined before. 

Comparing Eq.B.l to Eq.2, taking ~ 1, we see that 
a backwarming factor of Cbw = 2e -|- 1 is implicit in the 
IN05 calculation. Because Rrim oc \/Cbw, the IN05 subli- 
mation locations can be further than our estimates by up 
to a factor v'2e -|- 1 for any given grain type (to relate this 
to rim radii derived from interferometry, see Eq.5 and 6 
of Isella et al. 2006). Using Cbw(e) in Eq. A. 6 decreases 
differences with IN05 considerably (from 10% to 4% of 
Rings for 0.1/xm grains). However, as the comparison of 
results from such a Cbw said our Monte Carlo calculations 
has shown, this assumption generally overestimates the 
backwarming factor at the sublimation location. 

Midplane temperatures of a disc parametrized like an 
IN05 model with O.l/xm astronomical silicate are presented 
in Fig. B.2. Extrapolating the optically thin dust temper- 
ature curves (solid and dashed lines) to an intersection 
with the IN05 sublimation temperature curve for the same 
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density structure (dotted line), and taking into account 
the overestimated backwarming at the subhmation loca- 
tion, one sees how in MCMax, the inner rim is situated 
at 1.04AU and how it would be close to 1.16AU under 
the radiative transfer and dust sublimation assumptions 
of IN05. 

Monte Carlo radiative transfer was used by T07 to 
model discs parametrized like those of IN05, obtaining 
rim radii larger than the MCMax results in all cases, 
but smaller than IN05 for 1.2/im grains, and larger than 
IN05 for 0.1/xm. The TOT shift from MCMax seems to 
be clue mostly to their use of the IN05 dust sublima- 
tion law, i.e. a low Tgubb and a more basic treatment of 
dust sublimation and the associated numerical instabil- 
ities (A.Tannirkulam, personal communication; see also 
Appendix C). It is not clear whether the latter accounts 
for the changed order of rim radii for small and large grains 
between IN05 and TORUS (used by T07). 

In summary, compared to the treatment in the present 
work, the IN05 formalism leads to higher optically thin 
temperatures because of the assumptions implicit in the 
radiative transfer, and lower sublimation temperatures for 
pure olivine or astronomical silicate dust because of an 
inconsistent use of the dust density. These factors ex- 
plain their comparatively larger inner rim radii. The work 
of T07 adopted the IN05 sublimation law with its im- 
plications for the rim radii, and the treatment of subli- 
mation and condensation in TORUS is more basic than 
in MCMax, leading to further variations in rim location 
which we are presently unable to quantify. 

Appendix C: Numerical implementation 

In this appendix we discuss the numerical implementation 
of the sublimation physics described previously. We find 
that great care has to be taken in order to avoid instabil- 
ities or incorrect results. The iterative scheme described 
below is accompanied by a careful regridding of the spa- 
tial grid after each iteration to make sure that the optical 
depth through the disc is always sampled properly, allow- 
ing for accurate radiative transfer. Although the negative 
effects of a somewhat less optimized spatial grid are not 
always evident at first sight, we find that the influence 
on the temperature structure is significant and a proper 
spatial grid is of crucial importance for a proper imple- 
mentation of the sublimation physics. 



the inner region can already make the region in which this 
material was condensed optically thick, enhancing the ef- 
fects of backwarming and thereby lifting the temperature 
of the grains above the evaporation temperature again. At 
the same time, the region behind it is shielded and there 
dust can condense. This causes situations where the den- 
sity flips between two configurations which are both not 
the equilibrium solution. Thus, care has to be taken not 
to condense too nuich material in a single iterative step. 
However, the other extreme, condensing only so much ma- 
terial such that the change in optical depth is smaller than 
unity, requires a massive number of iterations which is un- 
feasible given the computation time required per iteration. 

Therefore, one must take care that the amount of ma- 
terial condensed or evaporated in a single iteration is not 
too small or too large. We have considered several differ- 
ent schemes for deciding how much dust to add or remove. 
The scheme we found to be the best tradeoff between sta- 
bility and speed decides how much dust to condense or 
evaporate according to how close the temperature of the 
dust species is to its evaporation temperature. In this way, 
it is possible to take large steps when the solution is far 
from equilibrium, while the size of the steps is automat- 
ically decreased when the solution is locally approached. 
The method is outlined below. 

The parameter we want to determine at every location 
in the disc is the gas fraction of eac;li spec;ies i. The 
gas fraction is defined as the fraction of the total available 
material for forming dust grains of species i that is in the 
gas phase. Thus 7^ = 1 means dust species i is totally 
evaporated, while 7i = means all the material is in the 
sohd phase. 

In order to determine 7,; at each loc;ation in the disc; 
we first compute the temperature structure for the initial 
guess of the dust density structure (see section C.2 be- 
low). Prom this we compute using Eq. (1) the partial gas 
pressure needed to counter the evaporation at each loca- 
tion in the disc. This gives the equilibrium gas fraction at 
this temperature, 7^,0- We also compute the sublimation 
temperature at each location, Tsubi and the temperature 
of the grains in the optically thin approximation, Tthin- 
The new gas fraction is then given by, 

7i,ncw = /u.7i,0 + (1 ^ /tu)7i,old, (C.l) 

where the weighting factor is taken as. 



C.l. Iterative method 

A straightforward iterative implementation of dust evapo- 
ration and recondensation does not always result in a sta- 
ble solution for the dust density structure. This is because 
of the high optical depths through the disc and the highly 
nonlinear radiative transfer effects resulting from it. In ad- 
dition, changing the dust density in one region of the disc 
affects the temperature structure at other places, making 
the problem highly non-local. For example, the conden- 
sation of a tiny fraction of the available dust material in 



T-T. 



subl 
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^subl 



for T < T,ubi 



for T > r,ubi 



(C.2) 

The parameter q determines how fast the condensa- 
tion of dust can take place and is adjusted such that the 
increase in optical depth up to the t ~ I surface as seen 
by direct stellar radiation is always less than 10% (i.e. 
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At < 0.1). The second term in the equation for for 
T < Tsubi ensures that in the region just outside the first 
few optical depths the condensation of material is not too 
fast. The temperature in these regions is quite low and if 
this term is not taken into account, the shielded regions 
will form so much dust that it will again start heating 
the region in front of it, causing another instability. By 
adding dust slowly also in these regions, the overall sta- 
bility of the method is increased significantly. Removal of 
dust grains when the temperature is too high is done much 
faster than recondensation. This ensures that during the 
iterative process only a tiny fraction of the dust grains has 
a temperature above the sublimation temperature. 

Together with the sublimation and recondensation of 
dust in the disc, we also determine the vertical density 
profile of the disc from hydrostatic equilibrium. 

C.2. Initial guess 

The success of most iterative schemes depend heavily on 
the initial guess. This is also true for the iterative scheme 
described above. Especially the estimate for the location 
of the start of the inner rim as a function of height above 
the midplane is of crucial importance for the success of 
the computation. We compute the initial guess for the dust 
density structure using analytic consideration as described 
below. 

First we take the temperature structure as obtained by 
Eq. (3) with an initial guess for Cbw (usually unity), ig- 
noring the thermal radiation and extinction of the stellar 
radiation by other dust grains. This temperature struc- 
ture can then be used to compute the gas fraction at each 
location in the disc, and from this we setup the first esti- 
mate for the density structure. We then do a full radiative 
transfer run through this density structure obtaining a 
temperature at each location. We use this to determine 
Cbw at the inner rim. Again we use Eq. (3) to determine 
the temperature at each location in the disc using the new 
value for Cbw This will result in a slightly different loca- 
tion of the inner rim. We iterate this procedure several 
times to get a proper value for Cbw and thereby a good 
estimate for the location of the inner rim. After this we 
proceed with the iterative scheme as described in the pre- 
vious section. 

C.3. Restrictions on the gas fraction gradients 

The location of the inner rim as a function of height in 
the disc is to a large extend determined by the local back- 
warming efficiency and the density of dust forming ma- 
terial. In the midplane of the disc the density is highest 
which causes the sublimation temperature to be highest 
and the rim radius to be small at this height. However, 
due to the heigh density the backwarming efficiency at the 
midplane is also highest, which causes the rim to retreat 
to larger radii. When this retreating effect of backwarming 
wins from the effect of a higher sublimation temperature. 



the rim radius at the midplane will be slightly larger than 
just above it causing a 'hole' in the disc. Inside this hole 
the dust is irradiated from almost all sides, causing the ef- 
fect of backwarming to increase even further and the rim 
to retreat even further. This is a runaway effect which, 
in the end, c;auses the entire disc to retreat to large radii. 
When this has happened, in front of this rim the tempera- 
tures are low and dust can start forming again. This cycle 
repeats and no stable solution is found. 

The above instability might be a physical effect which 
can emerge in a time-dependend treatment of the problem. 
However, we do not pursue this in this paper. Therefore, 
we restrict the gas fraction gradient to avoid triggering this 
instability. We restrict the gas fraction to be monotonously 
increasing with increasing height above the midplane. In 
this way we avoid the formation of a hole in the disc. 

Another problem which might arise from the iterative 
procedure explained above is that the destruction of dust 
proceeds faster in the region which is effectively back- 
warmed than in the optically thin region in front of it. 
This might cause the fraction of condensed material to 
decrease with increasing radius. This solution turns out 
to be unstable, and in addition we consider it to be an 
unlikely physical solution. Thus, to avoid this we restrict 
the gas fraction to be monotonously decreasing with in- 
creasing radius. 

Note that by restricting the gas fraction in the ways 
described above we force the rim to have the curved shape 
we see in the paper. However, other structures we have 
found have in all cases turned out to be unstable, thus we 
consider the curved rim a robust solution. 
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